\[\begin{align} \frac{dS}{dt} = & \left( \mu_s - \alpha_s - \alpha_{ss} S - \alpha_{ps}P \right) S + \beta_v V + \beta_p P \\ \frac{dP}{dt} = & \left( \alpha_s + \alpha_{ss} S \right) S - \left( \beta_p + \mu_p \right) P \\ \frac{dV}{dt} = & \left( \alpha_{ss} S + \alpha_{ps} P \right) S - \left( \beta_v + \mu_v \right) V \\ \frac{dO}{dt} = & \gamma_v V + \gamma_p P - \left( \delta_o + \mu_o \right) O \\ \frac{dC}{dt} = & \delta_o O - \mu_c C \\ \end{align}\]
The system of DEs is solved using LSODE. LSODE, the Livermore Solver for Ordinary Differential Equations, is a robust numerical solver known for its efficiency in handling stiff systems of differential equations. Stiff systems of ODEs are characterized by rapidly changing solutions, making them challenging to solve with standard methods. LSODE excels in tackling these complex systems. When dealing with a stiff system of ODEs, selecting the right numerical solver, such as LSODE, can significantly improve the accuracy and efficiency of the computations. To accurately simulate the behavior of the system, the initial values for the ODEs are extracted from the data recorded at the system’s initial time point.
Maximum Likelihood Estimation (MLE) is used to fit model parameters into data under the assumption that the error terms are Poisson distributed.
The data in question is a collection of crime reports and police dispatches from the Portland Police Bureau open data source. Note that we are only fitting the variables crime reports (\(C\)) and the police presence (\(O\)).
Let \(\widehat{O}\) and \(\widehat{C}\) be the estimated values of crime reports and police presence, respectively, from the model. Let \(o\) and \(c\) be the actual values of crime reports and police presence, respectively, from the data.
The goal of MLE is to determine the parameters for which the observed data have the highest joint probability. First, we write the parameters of the model as a vector
\[\vec{\theta} = \begin{bmatrix} \mu_s & \mu_p & \mu_v & \mu_o & \mu_c & \alpha_s & \alpha_{ss} & \alpha_{ps} & \beta_v & \beta_p & \gamma_v & \gamma_p & \delta_o \end{bmatrix}^T.\]
So, we let the function \(\vec{f}(t;\vec{\theta})\) be the solution to the system of differential equations. Since we only fit the variables \(O\) and \(C\) then
\[\vec{f}(t;\vec{\theta}) = \begin{bmatrix} \widehat{O} & \widehat{C} \end{bmatrix}^T.\]
Next, let \(X\) be a Poisson random variable with probability density function
\[P(X = x;\lambda) = \frac{\lambda^x}{x!} e^{-\lambda}\]
where \(\vec{\lambda} = \vec{f}(t;\vec{\theta})\) and \(\vec{x} = \begin{bmatrix} o & c \end{bmatrix}^T\). So, the likelihood function is written as
\[L(\lambda_1, \cdots, \lambda_n; x_1, \cdots, x_n) = \prod_{i=1}^{n} \frac{\lambda_i^{x_i}}{x_i!} e^{-\lambda_i}\]
where \(n\) is the length of the vector \(\vec{x}\). So, the negative log-likelihood is written as
\[\begin{align} -\ln{L(\lambda_1, \cdots, \lambda_n; x_1, \cdots, x_n)} = & -\ln{\left(\prod_{i=1}^{n} \frac{\lambda_i^{x_i}}{x_i!} e^{-\lambda_i}\right)} \\ = & - \sum_{i=1}^n \ln{\left(\frac{\lambda_i^{x_i}}{x_i!} e^{-\lambda_i}\right)} \\ -\ln{L(\lambda_1, \cdots, \lambda_n; x_1, \cdots, x_n)} = & - \sum_{i=1}^n \left( x_i \ln{\left( \lambda_i \right)} - \ln{\left( x_i ! \right)} - \lambda_i \right). \end{align}\]
We seek to find model parameter values \(\vec{\theta}\) to minimize/maximize the negative log-likelihood function written above. We used the Nelder-Mead optimization method to find the set of parameter values that minimizes the negative log-likelihood function.